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Several numerical relativity groups are using a modified ADM formulation for their simulations, 
which was developed by Nakamura et al (and widely cited as Baumgarte-Shapiro-Shibata-Nakamura 
system). This so-called BSSN formulation is shown to be more stable than the standard ADM 
formulation in many cases, and there have been many attempts to explain why this re-formulation 
has such an advantage. We try to explain the background mechanism of the BSSN equations by 
using eigenvalue analysis of constraint propagation equations. This analysis has been applied and 
has succeeded in explaining other systems in our series of works. We derive the full set of the 
constraint propagation equations, and study it in the flat background space-time. We carefully 
examine how the replacements and adjustments in the equations change the propagation structure 
of the constraints, i.e. whether violation of constraints (if it exists) will decay or propagate away. 
We conclude that the better stability of the BSSN system is obtained by their adjustments in the 
equations, and that the combination of the adjustments is in a good balance, i.e. a lack of their 
adjustments might fail to obtain the present stability. We further propose other adjustments to the 
equations, which may offer more stable features than the current BSSN equations. 



I. INTRODUCTION 

One of the current most important topics in the field 
of numerical relativity is to find a formulation of the 
Einstein equations which gives us stable and accurate 
longterm evolution. We all know that simulating space- 
time and matter based on general relativity is the es- 
sential research direction to go in the future, but we do 
not have a definite recipe for controlling numerical blow- 
ups. We concentrate our discussion on free evolution of 
the Einstein equations based on the 3-1-1 (space -I- time) 
decomposition of space-time, which requires solving the 
constraints only on the initial hypersurface and monitors 
the violation (error) of the calculation by checking con- 
straints during the evolution. 

Over the decades, the Arnowitt-Deser-Misner (ADM) 
[1] formulation has been treated as the default for nu- 
merical relativists. (More precisely, the version intro- 
duced by Smarr and York [2] was taken as the default, 
which we denote the standard ADM formulation here- 
after.) Although the ADM formulation mostly works for 
gravitational collapses or cosmological models in numer- 
ical treatments, it does not satisfy the requirement for 
longterm evolution e.g. the studies of gravitational wave 
sources. 

As we mentioned in our previous paper [3], we think 
we can classify the current efforts of formulating equa- 
tions for numerical relativity in the following three ways: 
(1) apply a modified ADM (BSSN) formulation [4, 5], (2) 
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apply a first-order hyperbolic formulation (see the refer- 
ences e.g. in [6, 7, 8]), or (3) apply an asymptotically 
constrained system [9, 10, 11, 12]. 

The first refers to using a modified ADM formulation, 
originally proposed by Nakamura in late 80s, and sub- 
sequently modified by Nakamura-Oohara and Shibata- 
Nakamura [4] . This introduces conformal decomposition 
of the ADM variables, a new variable for calculating Ricci 
curvature, and adjusts the equations of motion using 
constraints. The advantage of this formulation was re- 
introduced by Baumgarte and Shapiro [5] , and therefore 
this is often cited as the BSSN formulation, which we 
follow also. The BSSN equations are now widely used in 
the large scale numerical computations, including coales- 
cence of binary neutron stars [13] and binary black holes 
[14]. 

The second and third efforts use similar modifications 
such as introductions of new variables and/or adjust- 
ments of the equations, but differ in their purposes: to 
construct a hyperbolic formulation or to construct a for- 
mulation which constraints will decay or propagate away. 
The latter is intended to control numerical evolution so 
as the constrained manifold is its attractor. While the 
hyperbolic formulations have been extensively studied in 
this direction, we think the worrisome point in the discus- 
sion is the treatment of the non-principal part which is 
ignored in the hyperbolic formulation. As Kidder, Scheel 
and Teukolsky [8] reported recently, unless we reduce the 
effect of the non-principal part of the equations we may 
not get an advantages of hyperbolic formulation for nu- 
merical results [6, 15]. 

Through the series of studies [3, 6, 12, 16], we propose a 
systematic treatment for constructing a robust evolution 
system against perturbative error. We call it an asymp- 
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totically constrained (or asymptotically stable) system if 
the error decays itself. The idea is to adjust evolution 
equations using constraints (we term it an adjusted sys- 
tem), and to decide the coefficients (multipliers) by ana- 
lyzing constraint propagation equations. We propose to 
apply an eigenvalue analysis of the propagation equations 
of the constraints, especially in its Fourier components so 
as to include the non-principal part in the analysis. The 
characters of eigenvalues will be changed according to the 
adjustments to the original evolution equations. We con- 
jectured that the constraint violation which was occurred 
during the evolution will decay (if negative real eigenval- 
ues) or propagate away (if pure imaginary eigenvalues). 

This conjecture was affirmatively confirmed to ex- 
plain the numerical behaviors: wave propagation in the 
Maxwell equations [12], in the Ashtekar version of the 
Einstein equations [12], and in the ADM formulation (flat 
space-time background) [16]. The advantage of this con- 
struction scheme is that it can be applied to a formula- 
tion which is not a first-order hyperbolic form, such as 
to the ADM formulation [3, 16]. We think therefore our 
proposal is an alternative way to control/predict the vi- 
olation of constraints. (We believe that the idea of the 
constraint propagation analysis first appeared in Frittelli 
[17], where she made hyperbolicity classification for the 
standard ADM formulation). 

The purpose of this article is to apply this constraint 
propagation analysis to the BSSN formulation, and un- 
derstand how each improvement contributes to more sta- 
ble numerical evolution. Together with numerical com- 
parisons with the standard ADM case [18, 19], this topic 
has been studied by many groups with different ap- 
proaches. Using numerical test evolutios, Alcubierre et 
al. [20] found that the essential improvement is in the 
process of replacing terms by constraints, and that the 
eigenvalues of the BSSN evolution equations has fewer 
"zero eigenvalues" than those of ADM, and they conjec- 
tured that the instability can be caused by "zero eigenval- 
ues" that violate "gauge mode" . Miller [21] applied von 
Neumann's stability analysis to the plane wave propa- 
gation, and reported that BSSN has a wider range of 
parameters that give us stable evolution. These studies 
provide some supports regarding the advantage of BSSN, 
while it was also shown an example of an ill-posed solu- 
tion in BSSN (as wefl in ADM) [22]. (Inspired by BSSN's 
conformal decomposition, several related hyperbolic for- 
mulations have also been proposed [23, 24, 25].) 

We think our analysis will offer a new vantage point on 
the topic, and contribute an alternative understanding of 
its background. Consequently, we propose more effective 
improvement of the BSSN system which has not yet been 
tried in numerical simulations. 

The construction of this paper is as follows. We review 
the BSSN system in §11, and there also we discuss where 
the adjustments are applied. In §111 we apply our con- 
straint propagation analysis to show how each improve- 
ment works in the BSSN equations, and in §IV we ex- 
tend our study to seek a better formulation which might 



be obtained by small steps. We only consider the vac- 
uum space-time thoughout the article, but the inclusion 
of matter is straightforward. 



II. BSSN EQUATIONS AND THEIR 
CONSTRAINT PROPAGATION EQUATIONS 

A. BSSN equations 

We start presenting the standard ADM formulation, 
which expresses the space-time with a pair of 3-metric 
7ij and extrinsic curvature Kij . The evolution equations 
become 

d^j.j = -2aK,j+D,f3j+Djl3,, (2.1) 

+{Dil3'')Kk, + {D.j0'')Kk^ + P^DkK.j (2.2) 

where a, (3i are the lapse and shift function and Di is the 
covariant derivative on 3-space. The symbol means 
the time derivative defined by these equations, and we 
distinguish them from those of the BSSN equations , 
which will be defined in (2.15)-(2.19). The associated 
constraints arc the Hamiltonian constraint Ti and the 
momentum constraints Aii'. 

j^ADM ^ rADM ^j^2 _J^,,J^^J^ (2.3) 

^ADM ^ DjKU-D,K. (2.4) 

The widely used notation[4, 5] is to introduce the vari- 
ables {^p,^ij,K,Aij,T^) instead of {-fij,Kij), where 





= (l/12)log(det7i,). 


(2.5) 




= e"'''^7ij. 


(2.6) 


K 




(2.7) 


Aij 


= e-^^{Ki, - {l/iHjK), 


(2.8) 


pi 




(2.9) 



The new variable F* was introduced in order to calculate 
Ricci curvature more accurately. F' also contributes to 
make the system re-produce wave cqiiations in its linear 
limit. In the BSSN formulation, Ricci curvature is not 
calculated as 

Rfj^'^' = i9feFy - diT^j + F-^Fj\ - T^jTfi, (2.10) 

but 

RBSSN ^ RP^+R^^^ (2.11) 
Rf. = -2b,bj'f-2y,jb''bk^ 

+A{b,^){b,ip) - Aj,jib''ip){bk^), (2.12) 

Rij = -(l/2)7"^9,afe%-+7fc(,%f'=+f'=f(,,)fe 

+2f ™f f(,f,.),„ + ¥"TL^kij, (2.13) 
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where Di is covariant derivative associated with "/ij. 
These are weakly equivalent, but Rf^^^ does have wave 
operator apparently in the flat background limit, so that 
we can expect more natural wave propagation behavior. 

Additionally, BSSN requires us to impose the confor- 
mal factor as 

7(:= det%) = 1, (2.14) 
during the evolution. This is a kind of definition, but can 



also be thought of as a constraint. We will return to this 

point shortly. 

BSSN's improvements are not only the introductions 
of new variables, but also the replacement of terms in the 
evolution equations using the constraints. The purpose 
of this article is to understand and to identify which im- 
provement works for the stability. Before doing that we 
first show the standard set of the BSSN evolution equa- 
tions: 



d^^ = -{l/6)aK + {l/6)/3\d,^) + (a,/3''), (2.15) 

d^jij = -2aA,, + 7,fe(a,/3'=) + jjkidrP") - (2/3)7^(9^/3'=) + /3\dk%), (2.16) 

afif = ~D'D^a + aA,jA'^ +{l/3)aK^ + (3'{^^K), (2.17) 
dfAij = -e-^'^iDiD.af^ + e-^'^a(i?,^^^^)^^ + aKAij - 2aiifei^- + (9^/3'=)^, + (a,/3'=)ifei 

-{2/3){dk0'')Aij + p\dkAij), (2.18) 
9f P = -2{dja)A'^ + 2a(f jfei'^^' - {2/3)f^{djK) + 6A'^idj<fi)) - dj{(3\dkl'^) - j'^idkp') 

-l''\dk/3^) + {2/3)f^{dkf3'')). (2.19) 



We next summarize the constraints in this system. The 
normal Hamiltonian and momentum constraints (the 
"kinematic" constraints) are naturally written as 



n 



BSSN 



R 



BSSN 



K^-KijK'\ 



(2.20) 
(2.21) 



where we use Ricci scalar defined by (2.11). Addition- 
ally, we regard the following three as the constraints (the 
"algebraic" constraints): 



(2.22) 



A 



5 = 7-1, 



(2.23) 
(2.24) 



where the first two are from the algebraic definition of 
the variables (2.8) and (2.9), and the (2.24) is from the 
requirement of (2.14). Hereafter we write 'H^^^^ and 
j^BSSN gjjjjpiy as TL and M. respectively. 

Taking careful accomit of these constraints, (2.20) and 
(2.21) can be expressed directly as 



n = e'^'fR - 8e-^'^DWjip - Se-'^'^{& ^){bjip) + (2/3)/^^ - AijA^^ - {2/2,)AK, (2.25) 
Mi = &Ah{Dj^)-2A{D,ip)-{2/?>){D,K)+^''^{b^Ak^). (2.26) 

In summary, the fundamental dynamical variables in BSSN are {(p,%j, K,Aij,r^), total 17. The gauge quantities 
are (a, /3*) which is 4, and the constraints are {H, A4i, G^, A, S), i.e. 9 components. As a result, 4 (2 by 2) components 
are left which correspond to two gravitational polarization modes. 



B. Adjustments in evolution equations 

Next, we show the BSSN evolution equations (2.15)-(2.19) again, identifying where the terms are replaced using 

the constraints, (2.20)-(2.24). 

By a straightforward calculation, we get: 

9f = dfip+il/6)aA~il/12)r\d,S)f3', (2.27) 
= - {2/3)a%A+{l/3)r\dkS)p''%, (2.28) 
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d^K = dfK-{2/3)aKA-an + ae~'^'^{D^g^), (2.29) 
5f ii,- = d^A, + {{l/3)a%K - {2/3)aAj)A+{{l/2)ae-^^{dk7ij) - {l/6)ae-^^%r\dkS))g'' 

+ae-4'^7M.(9,)a') - {l/3)ae-*'^%{dkg'), (2.30) 
r = dfr + ( - (2/3)(aja)7^'* - i2/3)a{djf') - {l/3)aj''T\d,S) + Aaf^d,ip))A - i2/3)af'{djA) 

+{5/6)p'^r^'^idkS)id,S) + il/2)p''r\dkfnid,S) + {l/3)p'^r\djjn{dkS). (2.31) 



where denotes the part of no replacements, i.e. the 
terms only use the standard ADM evolution equations in 
its time derivatives. 

From (2.27)-(2.31), we understand that all the BSSN 
evolution equations are adjusted using constraints. This 
fact will give us the importance of the scaling constraint 
<S = and the tracefree operation A = during the 
evolution. 

As we have pointed out in the case of adjusted ADM 
systems [3, 16], certain combinations of adjustments (re- 
placements) in the evolution equations change the eigen- 
values of constraint propagation equations drastically. 
For example, all negative eigenvalues can be negative 
real by applying Detweiler's adjustment [26] or its sim- 
plified version. One common fact we found is that such a 
case has an adjustment which breaks time reversal par- 
ity of the original equation. That is, with a change of 
time integration direction dt — > —dt, an adjusted term 
might become effective if it breaks time reversal symme- 
try. (This time asymmtric feature was first imprementcd 
as a "lambda-system" in [9].) Unfortunately, for the case 
of the BSSN equations, (2.27)-(2.31), all the above ad- 
justments keep the time reversal symmetry. So that we 
can not expect direct decays of constraint violation in the 
present form. We will give the details on this point later. 

III. CONSTRAINT PROPAGATION ANALYSIS 
IN FLAT SPACE-TIME 

A. Procedures 

We start this section overviewing the procedures and 

our goals. In our scries of previous works [3, 12, 16], we 
have concluded that eigenvalue analysis of the constraint 
propagation equations are quite useful for explaining or 
predicting how the constraint violation grows. 

Suppose we have a set of dynamical variables it"(x*, t), 
and their evolution equations 

atU« = /(M«,9i««,---), (3.1) 

and the (first class) constraints 

C"(u^aiU^•••) «0. (3.2) 

For monitoring the violation of constraints, we propose 
to investigate the evolution equations of C" (constraint 
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propagation), 

dtC"=g{C",diC",---). (3.3) 

(We do not mean to integrate (3.3) numerically, but 
rather to evaluate them analytically in advance.) In order 
to analyze the contributions of all RHS terms in (3.3), we 
propose to reduce (3.3) in ordinary differential equations 
by Fourier transformation, 

dtC"" = 3(6") = M^/sCf^, (3.4) 

where C{x,t)P = JC{k,t)P cxp(ik ■ x)d^k, and then to 
analyze the set of eigenvalues, say Aq,, of the coefficient 
matrix, M"/}, in (3.4). We call As and M"^ the con- 
straint amplification factors (CAFs) of (3.3) and con- 
straint propagation matrix, respectively. Our guidelines 
to have 'better stability' are that 

(A) If the CAFs have a negative real-part (the con- 
straints are forced to be diminished), then we see 
more stable evolution than a system which has a 
positive real-part. 

(B) If the CAFs have a non-zero imaginary-part (the 
constraints are propagating away), then we see 
more stable evolution than a system which has zero 
CAFs. 

We found heuristically that the system becomes more 
stable when more As satisfy the above criteria [6, 12]. 
We note that these guidelines are confirmed numerically 
for wave propagation in the Maxwell system and in the 
Ashtekar version of the Einstein system [12], and also 
for error propagation in Minkowskii space-time using ad- 
justed ADM systems [16]. Supporting theorems for above 
(A) was recently discussed [31]. 

The above features of the constraint propagation, 
(3.3), will differ when we modify the original evolution 
equations. Suppose we add (adjust) the evolution equa- 
tions using constraints 

dtu'' = /(««, a^n", ■ • •) + F(C«, a,C-, • • •), (3.5) 

then (3.3) will also be modified as 

atC« = g{C", diC, ■■■)+ G((7«, diC, ■ ■ •). (3.6) 

Therefore, the problem is how to adjust the evolution 

equations so that their constraint propagation satisfies 
the above criteria as much as possible. 
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B. BSSN constraint propagation equations 

Our purpose in this section is to apply the above pro- 
cedure to the BSSN system. The set of the constraint 
propagation equations, Mi, G^, A, 5)^, turns to be 

quite long and not elegant (is not a first-order hyperbolic 
and includes many non- linear terms), and we put them 



in Appendix. In order to understand the fundamental 
structure, we hereby show an analysis on the flat space- 
time background. 

For the flat background metric g^i^ = rj^^, the first or- 
der perturbation equations of (2.27)-(2.31) can be written 
as 



dPcP = -(l/6f k + (1/6)(k^ - if U, (3.7) 

dt^%j = -i^^Aij - {2/3){k^ - l)SiPA, (3.8) 

= -{djd/^h) + {KKl-l)d/^^^ -{KK2-lf^, (3.9) 

a/D^. = ^^~)^^ -(^\D,D,af^ + {kai - 1)4(^(9,/^^'=) - {1/3){ka2 - mj{dk^'^% (3.10) 

= -(4/3)(5/i)i^) - (2/3)(Kfi - l){dPA) + 2(Kf2 - ^f^Mi, (3.11) 

where we introduced parameters ks, all k = reproduce no adjustment case from the standard ADM equations, and 
all K = 1 correspond to the BSSN equations. We express them as 

(3.12) 

Constraint propagation equations at the first order in the flat space-time, then, become: 

9/1)^ = (k;^ - (2/3)/tpi - (4/3)k^ + 2) djdj'^^U + 2(Kf2 - 1)(6>/^)Mj), (3.13) 
dPMi = {-{2/3)kki + (l/2)«Ai - (1/3)ka2 + (1/2)) didj^^^^ 

+{l/2)KAidjdj^^^' + {{2/3)kk2 - (1/2)) 9/%, (3.14) 

a/i^* = 2Kf2(ilV(,: + (-(2/3)Kn - (l/3)/«;^)(a/iU), (3.15) 

= -2K^(il4, (3.16) 

dPA = (k^i-«A2)(5/i^^). (3.17) 

We win discuss CAFs of (3.13)-(3.17). 



C. Effect of adjustments 

We check CAFs of the BSSN equations in detail. The 
list of examples is shown also in Table I. Hereafter we 
let k"^ = k'i -\- k'^, + fc^ for Fourier wave numbers. 

1. The no-adjustment case, Kadj =(all zeros). This is 
the starting point of the discussion. In this case, 

CAFs = (0(x7),±\/^), 

i.e., (0 (x 7), ibpure imaginary (1 pair)). In the 
standard ADM formulation, which uses {jij,Kij), 
CAFs arc (0, 0, ±Pure Imaginary) [16]. Therefore 
if we do not apply adjustments in the BSSN equa- 
tions the constraint propagation structure is quite 
similar to that of the standard ADM. 



i.e., (0 (x 3), ibPure Imaginary (3 pairs)). The 
number of pure imaginary CAFs is increased over 
that of No.l, and we conclude this is the advantage 
of adjustments used in the BSSN equations. 

3. No 5-adjustment case. All the numerical experi- 
ments so far apply the scaling condition S for the 
conformal factor ip. The tS-originated terms appear 
many places in the BSSN equations (2.15)-(2.19), 
so that we suspect non-zero 5 is a kind of source 
of the constraint violation. However, since all S- 
originated terms do not appear in the flat space- 
time background analysis, [no adjusted terms in 
(3.7)-(3.11)], our analysis is independent of the <S- 
constraint. (Remark that we do not deny the effect 
of <S-adjustment in other situation.) 



2. For the BSSN equations, Kadj =(all Is), 



CAFs = (0 (x3), ±\/^ (3 pairs)). 



4. No .A-adjustment case. The trace (or traceout) 

condition for the variables is also considered nec- 
essary (e.g. [27]). This can be checked with 



6 



Kadj = (k, k, 1, 1, 1, 1, K, 1), and we get 
CAFs = (0(x3),±V^(3 pairs)), 

independent of k. Therefore the effect of A- 
adjustment is unimportant according to this anal- 
ysis, i.e. on flat space-time background. (Remark 
that we do not deny the effect of ^-adjustment in 
other situation.) 

5. No ^*-adjustmcnt case. Tlio introduction of F* is 
the key in the BSSN system. If we do not apply ad- 
justments by Q^, {Kadj = (1,1,0,1,0,0,1,1)) then 
we get 

CAFs = {0{x7),±^/^), 

which is the same with No.l. That is, adjustments 
due to terms are effective to make a progress 
from ADM. 

6. No A^i-adjustment case. This can be checked with 

f^adj = (1. 1, 1, 1, 1, 1, 1, k), and we get 

CAFs = (0,±V-Kfc2 (2 pairs), 

±y/-k^-l + 4K+\l-4K\)/6, 

±V'-A;2(-1 + 4:K-\1- 4k|)/6). 

If K = 0, then (0(x7),±v/fcV3), which is 
(0(x7), ibreal value). Interestingly, these real val- 
ues indicate the existence of the error growing mode 
together with the decaying mode. Alcubierre et al. 
[20] found that the adjustment due to the momen- 
tum constraint is crucial for obtaining stability. We 
think that they picked up this error growing mode. 
Fortunately at the BSSN limit {k = 1), this error 
growing mode disappears and turns into a propa- 
gation mode. 

7. No H-adjustment case. The set Kadj = 
(1,1,1,K,1,1,1,1) gives 

CAFs = (0 (x3), ±\/^ (3 pairs)), 

independently to k. Therefore the effect of H- 
adjustment is unimportant according to this anal- 
ysis, i.e. on flat space-time background. (Re- 
mark again that we do not deny the effect of H- 
adjustment in other situation.) 

These tests are on the effects of adjustments. We will 
consider whether much better adjustments are possible 
in the next section. 

We list the above results in Table I. (Table I includes 

a column of diagonalizability of constraint propagation 
matrix M, of which importance was pointed out in [31]. 
) The most characteristic points of the above are No. 5 
and No. 6 that denote the contribution of the momentum 
constraint adjustment and the importance of the new 



variable F*. It is quite interesting that the unadjusted 
BSSN equations (case 1) docs not have apparent advan- 
tages from the ADM system. As we showed in the case 
5 and 6, if we missed a particular adjustment, then the 
expected stability behavior occationally gets worse than 
the starting ADM system. Therefore we conclude that 
the better stability of the BSSN formulation is obtained 
by their adjustments in the equations, and the combina- 
tion of the adjustments is in a good balance. That is, a 
lack of their adjustments might fail to obtain the stability 
of their system. 



IV. PROPOSALS OF IMPROVED BSSN 
SYSTEMS 

In this section, we consider the possibility whether 
we can obtain a system which has much better prop- 
erties; whether more pure imaginary CAFs or negative 
real CAFs. 



A. Heuristic examples 

(A) A system which has 8 pure imaginary CAFs: 

One direction is to seek a set of equations which make 
fewer zero CAFs than the standard BSSN case (No. 2 in 
the previous section). Using the same set of adjustments 
in (3.7)-(3.11), CAFs are written in general 

CAFs = (o,±y^ 

'«AiKf2 (2 pairs), 
ibcomplicated expression, 
ibcomplicated expression^ . 

The terms in the first line certainly give four pure imag- 
inary CAFs (two positive and negative real pairs) if 
K^iKp2 > (< 0). Keeping this in mind, by choosing 
i^adj = (1, 1, 1, 1, 1, K, 1, 1), we find 

CAFs = (0,±\/QI^(2 pairs). 

±v/-A-^(2 + K + |/v-4|)/6, 
±^/-k^{2 + K-\K-4\)/6,y 

Therefore the adjustment Kadj = (1, 1, 1, 1, 1, 4, 1, 1) gives 

CAFs = (o,±^/^ (4 pairs)), 

which is one step advanced from BSSN's according our 
guidelines. 

We note that such a system can be obtained in many 
ways, e.g. Kadj = (0,0.1,0,2,1,0,1/2) also gives four 
pairs of pure imaginary CAFs. 

(B) A system which has negative real CAF: 

One criterion to obtain a decaying constraint mode (i.e. 
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an asymptotically constrained system) is to adjust an 
evolution equation as it breaks time reversal symmetry [3, 
16]. For example, we consider an additional adjustment 
to the BSSN equation as 



(4.3) 



(4.1) 



which is a similar adjustment of the simplified Detweiler- 
type [26] that was discussed in [3]. The first order con- 
straint propagation equations on the flat background 
space-time become 

dPMi = {l/Q)dl^M+{l/2)djdp-^\ 
dPA = -{d^d.^^hf^ + <^^^Rff^''f^, 

= -i')A + iKSD^'% 

where we wrote only additional terms to (3.13)-(3.17). 
The CAFs become 

CAFs = (0(x2), 

±v/^(3 pairs), {3/2)k^KsD), 

in which the last one becomes negative real if ksd < 0. 

(C) Combination of above (A) and (B) 
Naturally we next consider both adjustments: 



(4.2) 



J 



where the second one produces the 8 pure imaginary 
CAFs. The additional terms in the constraint propa- 
gation equations (3.13)-(3.17) are 



= djd/^U- {3/2)KsDdjdj^^^, 

dPMi = {i/Q)dl^M + {i/2)djdl-^^' 

5/1^^ = -dPA+ {l/2)KsDdl^^ + i^^Mi, 
a/^U = -3«85fc(^^^ 

We then obtain 



CAFs = (o,±V^(3 pairs). 



{2,/A)k'KSD ± ^F(-K8 + (9/16)A;24^)) 



which reproduces case (A) when ksd = 0, ks = 1, and 
case (B) when Kg = 0. These CAFs can become (0, pure 
imaginary (3 pairs), complex numbers with a negative 
real part (1 pair)), with an appropriate combination of 
Kg and KSD- 



B. Possible adjustments 

In order to break time reversal symmetry of the evolution equations [3, 9, 16], the possible simple adjustments are 
(1) to add W, <S or Q'^ terms to the equations of dtcf), 9t7ij, or dtV^, or (2) to add M.i or A terms to dtK or dtAij. We 
write them generally, including the above proposal (B), as 



dt(t> = df(j) + K^uoi^ + HQ'^^kG^-> 
dtlij = dflij + Ksd a^H + K^gi a^DkO'' + K^g2 Oi^k(iDj)G'^ + k^si a^ijS + aDiDjS, 
dtK = dfK + KKMaf''{D,Mk), 
dtA^j = dfAij + KAMI ctlij {D^Mk) + KAM2 a{D(^iMj)) + kaai ajijA + kaa2 aDiDjA, 
= dfr + Kf^ aD'H + Kf aQ' + Kfg^ aDWjQ' + Kf,g^ aD'DjQ^ , 

where ks are possible multipliers (all k = reduce the syste m the standard BSSN evolution equations). 



(4.4) 
(4.5) 
(4.6) 
(4.7) 
(4.8) 



We show the effects of each terms in Table II. The 
CAFs in the table are on the flat space background. We 
see several terms produce negative real-part in CAFs, 
which might improve the stability than the previous sys- 
tem. (Table II includes again a column of diagonalizabil- 
ity of constraint propagation matrix M. Diagonalizablc 
ones are expected to reflect the predictions from eigen- 
value analysis. That is, the eigenvalue analysis with diag- 



onalizablc ones definitely avoid the diverging possibility 
in constraint propagation when it includes degenerated 
CAFs. Sex; [31]). For the readers convenience, we list up 
several best candidates here. 

(D) A system which has 7 negative CAFs 
Simply adding D^^^Aij) term to dtAij equation, say 

dtAij = df^^^Aij + KAM2a{D^iMj) ) (4.9) 
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with K,AM2 > 0, CAFs on the flat background are 7 
negative real CAFs. 

(E) A system which has 6 negative CAFs 
The below two adjustments will make 6 negative real 
CAFs, while they also produce one positive real CAF (a 
constraint violating mode). The efl'ectiveness is not clear 
at this moment, but we think they are worth to be tested 
in numerical experiments. 

(El) With K^g2 < 0, 

daij = df^^^^ij + K^g2 a^kiiDj) ■ (4.10) 
(E2) With Kfg^ < 0. 

d^r = d^^^^r + Kfg^ a&bjg\ (4.ii) 

V. CONCLUDING REMARKS 

Applying the constraint propagation analysis, we 
tried to understand why and how the so-called BSSN 
(Baumgartc-Shapiro-Shibata-Nakamura) re- formulation 
works better than the standard ADM equations in gen- 
eral relativistic numerical simulations. Our strategy was 
to evaluate eigenvalues of the constraint propagation 
equations in their Fourier modes, which method suc- 
ceeded to explain the stability properties in many other 
systems in our series of works. 

We have studied step-by-step where the replacements 
in the equations affect and/or newly added constraints 
work, by checking whether the error of constraints (if it 
exists) will decay or propagate away. Alcubierre et al [20] 
pointed out the importance of the replacement (adjust- 
ment) to the evolution equation using the momentum 
constraint, and our analysis clearly explains why they 
concluded this is the key. Not only this adjustment, we 
found, but also other adjustments and other introduc- 
tions of new constraints also contribute to making the 
evolution system more stable. We found that if we missed 
a particular adjustment, then the expected stability be- 
havior occationally gets worse than the ADM system. We 



further propose other adjustments of the set of equations 
which may have better features for numerical treatments. 

The discussion in this article was only in the flat back- 
ground space-time, and may not be applicable directly to 
the general numerical simulations. However, we rather 
believe that the general fundamental aspects of con- 
straint propagation analysis are already revealed in this 
article. This is because, for the ADM and its adjusted 
cases, we found that the better formulations in the flat 
background are also better in the Schwarzschild space- 
time, while there are differences on the effective adjusting 
multipliers or the effective coordinate ranges[3, 16]. 

We have not shown any numerical tests here. However, 
recently, the proposal (B) in §IV was examined numeri- 
cally using linear wave initial data and confirmed to be 
effective for controlling the violation of the Hamiltonian 
constraint with our predicted multiplier signature [28]. 
The systematic numerical comparisons between different 
formulations are underway [29], and we expect to have a 
chance to report them in near future. We are also try- 
ing to explain the stability of Laguna-Shoemaker's impli- 
mented BSSN system [30] using the constraint propaga- 
tion analysis. 

There may not be the almighty formulation for any 
models in numerical relativity, but we believe our guide- 
lines to find a better formulation in a systematic way 
will contribute a progress of this field. We hope the pre- 
dictions in this paper will help the community to make 
further improvements. 
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APPENDIX A: FULL SET OF BSSN CONSTRAINT PROPAGATION EQUATIONS 

The constraint propagation equations of the BSSN system can be written as fohows. 
dtH = ({2/i)aK + {2/i)aA + l3''dk)n+ ( - Ae-^"^ a{dkip)^''^ - 2e-^'^{dka)f''^Mj 

+ (2ae-4'^7-i7"=(a,^)^afc + {l/2)ae-^^r\diAW''dk + {l/2)e-^^r\Qia)l"' Ad, 

+(l/2)e-^^7-'/3'(9,-5iy')9fc + {'i/A)e-^'^r''P'l'\diS){djS)dk - {^/A)e-^^r^(3'{da^''){djS)dk 
+(l/3)e-4'^7-'7^'(9,/3'=)apafc 

-{^/12)e-^^r^^^\dkP'){diS)dj + {l/^)e-^'^r\dkl'^){djfi'')di - {l/6)e-^'^r^^"'\dkdip')dm)s 

+ ({A/%)aKA-{S/9)aK^ + {A/i)ae-^^ {d^dj^)f' + {8/?,)ae-^'P{dk'p){dn"') 

+2e-^^{dia)^'''dk + e-^^^'\didka))A, (Al) 
dtM, = ( - (l/3)(a,a) + (l/Q)d^n + aKM, + (ae-^^^^^^{dkv){dj^rm) - (l/2)ae-4^fS7''(5,7™,) 

+ (l/2)ae-4^7"'^(afe9,7™) + {l/2)ae-^^r\d,S){d,S) - {l/A)ae-^'^{d,%i){d,^''') + ae-'"^^'''^{dkV>)ljidm 
+ae-^^(a,^)9, - {l/2)ae-^'^rTlt'lJ^^m + ae-^^f "'=f„-fea„ + (l/2)ae-^'^j"'j,,dkdi 
+{^/2)e-^^j"'\dj^im){dka) + {l/2)e-^^{dja)di + {l/2)e-'"^r''^j,,{dka)dm)g' 

+ ( - i^C^fca) + {l/9){dia)K + {A/9)a{diK) + {l/9)aKdi - ai^Qfc) A (A2) 

dtG' = 2af^Mj + ( - {l/2)p'' f'r\diS)dk - {l/2)(}''f^{dk^mn)7'^'r'di + {l/2)(3''f'r'didk 

-{l/2){d^(3')^""'j-^dk + {l/3){dip')f''r'dk)s + ( + Aaf^Dj^) - af'd, - {dkaW'^^A, (A3) 
dtS=+l3''idkS)-2ajA, (A4) 
dtA= (aK + p''dk)A. (A5) 
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The flat background linear order equations, (3.13)-(3.17), were obtained from these expression. 



No. 




Constraints (number of components) 


diag? 


CAFs 


in text. 




n{i) 


Mi (3) g' (3) A (1) 


5(1) 




in Minkowskii background 


0. 


standard ADM 


use 


use 


- 


yes 


(0,0,9,9) 


1. 


BSSN no adjustment 


use 


use use use 


use 


yes 


(0, 0, 0, 0, 0, 0, 0, 9, 9) 


2. 


the BSSN 


use+adj 


use+adj use+adj use+adj 


use+adj 


no 


(0, 0, 0, 9, 9, 9, 9, 9, 9) 


3. 


no S adjustment 


use+adj 


use+adj use+adj use+adj 


use 


no 


no difference in flat background 


4. 


no A adjustment 


usc+adj 


use+adj use+adj use 


use+adj 


no 


(0,0,0,9,9,9,9, 9, 9) 


5. 


no Q'' adjustment 


use+adj 


use+adj use use+adj 


use+adj 


no 


(0,0,0,0,0,0,0,9,9) 


6. 


no Mi adjustment 


use+adj 


use use+adj use+adj 


use+adj 


no 


(0,0, 0,0,0,0, 0, »,») 


7. 


no Ti. adjustment 


use 


use+adj use+adj use+adj 


use+adj 


no 


(0, 0, 0, 9, 9, 9, 9, 9, 9) 



TABLE I: Summary of §111 C: contributions of adjustments terms and effects of introductions of new constraints in the BSSN 
system. The center column indicates whether each constraints are taken as a component of constraints in each constraint 
propagation analysis ('use'), and whether each adjustments are on ('adj'). The column 'diag?' indicates diagonalizability of 
the constraint propagation matrix. The right column shows CAFs, where 9 and 3? means pure imaginary and real eigenvalue, 
respectively. No.O (standard ADM) is shown in [16]. 



adjustment 


dt<f) 


K^n oiH 


dt<t> 






KsD a'jij'H 








K^g2 a^k(iDj)Q'^ 


dtlij 


Kxisi a^ijS 


dtlij 


K752 aDiDjS 


dtK 


kkm a^^''{DjMk) 


dtAij 


KAMI ajij(D''Mk) 


dtAij 


KAM2 a{D(^iMj)) 


dtAij 


KAAi a^ijA 


dtAij 


KAA2 aDiDjA 


dtr 


aD^H 


dtr 


Kf 61 


dtV' 


K^g^a&Djg' 


dtV' 


Kf^g^a&Djg^ 



CAFs 



(0,0,±V^(*3),8«:^wfc2) 

(0, 0, ±V— fc^(*2), long expressions) 

(0, 0, ±V^i*3), (3/2)/tsDfc^) 

(0, 0, ±V— fc'^(*2), long expressions) 

(0,0, (l/4)fe2«;-g2 ± y'A;2(-l+fe2K.^g2/16)(*2), 

long expressions) 
(0,0,±\/^(*3),3«:^si) 

(0,0,0,±y^(*2) 



dias? 



effect of the adjustment 



K.4>H < makes 1 Neg. 

ii4,g < makes 2 Neg. 1 Pos. 

i^SD < makes 1 Neg. Case (B) 

Kjgi > makes 1 Neg. 

^762 < makes 6 Neg. 1 Pos. Case (El) 

lijSi < makes 1 Neg. 
ii^S2 > makes 1 Neg. 

i^KM < makes 2 Neg. 

i^AMi > makes 1 Neg. 

iiAM2 > makes 7 Neg Case (D) 

i^AAi < makes 1 Neg. 

KAA2 > makes 1 Neg. 

Kp^ > makes 1 Neg. 

i^tgi < makes 6 Neg. 1 Pos. Case (E2) 

Kfgj > makes 2 Neg. 1 Pos. 

Kpg3 > makes 2 Neg. 1 Pos. 



)) 



J, u, u, inv 

(l/3)KxMfc2 ± (l/3)^A;2(-9 + fc2K2,^ 
(0,0,±v^(*3),-KAAiifc^) 

(0,0, -fc^KA.VI2/4 ± + fc2KAM2/16)(*2) 

long expressions) 
(0,0,±^/^(*3),3/tA.4l) 

(0,0,±y^(*3),-«A.A2fc^) 
(0, 0, ±V^{*S), - KAA2k^) 

(0,0, {l/2)Kpg, ± ^/-F + K?g^(*2 



long cxp.) 

(0,0,-(l/2)Kfg2 ± yiF+4^(*2) , long exp. 
(0,0, -(l/2)Kfe3 ± J-k''+4g^i*2) , long exp. 



no 

yes 
yes 
yes 

yes 

no 
no 

no 

yes 
yes 

yes 
yes 
no 
yes 

yes 

yes 



TABLE 11: Possible adjustements which make a real-part CAFs negative (§IV B). The column of adjustments are nonzero 
multipliers in terms of (4.4)-(4.8), which all violate time reversal symmetry of the equation. The column 'diag?' indicates 
diagonalizability of the constraint propagation matrix. Neg./Pos. means negative/positive respectively. 



